/*
 * Analyze_b.h
 *
 *  Created on: Mar 1, 2012
 *      Author: marchi
 */

#ifndef ANALYZE_B_H_
#define ANALYZE_B_H_


#include "Metric.h"
#include "CellSymmetry.h"
#include "Voronoi.h"
#include "gmx_voronoi.h"
#include "Field.h"
#include "Parameters.h"
#include "gmx_voroelec.h"
#include "Gridn_SpecTemplates.hpp"
#include "Grid_SpecTemplates.hpp"
#include "Communicator.hpp"

#include "DiffCoeffs_Spec.hpp"
#include "LocalEps.h"
#include <cmath>

using std::vector;
using namespace gmx_Voronoi;
using namespace gmx_voroelec;
const bool TEST=true;
static int nframes=0,nstart,freq;
static int myframes=0;
static CellSymmetry * symmetry=NULL;
gmx_ana_nbsearch_t  * Mynb=NULL;
static bool bHydrogen=FALSE;
static double dtt;
static double Mytime=0.0;

static int Analyze(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{

    /* Here, you can do whatever analysis your program requires for a frame. */
	GridEps & MyEps=*MyEps0;
    Voronoi & vor=*(static_cast<Voronoi *>(data));

    MyDipole & Dip=*Dips;
	if(!symmetry)
		symmetry=new CellSymmetry(top,fr,pbc);

	symmetry->Center(nindex,cindex, natoms);
    Metric Met;
    Met(fr->box);
    Grid<1>::set(fr->box);
	Gridn<DIM>::set(fr->box);

    Dip.setCoord(Met, fr->x);

    P_t->Density(Dip);
    P_t->Filter();

    Rho_t->Density(Dip);
   (*E_t)(*Rho_t);

    if(!MyEps(*P_t, *E_t,false)) exit(1);
    if(nframes == 0 && Parameters::Input::freq > 0){
    	dtt=static_cast<double> (fr->time)/static_cast<double> (fr->step);
    	nstart=fr->step;
    	freq=rint(Parameters::Input::freq/dtt/Parallel::comm->Get_Size());
    	Parallel::comm->Broadcast(&freq,1);
    }
    if(nframes > 0 && Parameters::Input::freq > 0) {
    	int actual=fr->step-nstart;
    	Parallel::comm->Broadcast(&actual,1);

    	if(!(actual % freq)){
    		myframes++;
    		MyEps.setMetric(fr->box);
    		MyEps.setxRef(Dvect(0.0,0.0,0.0));
    		fout << MyEps ;
    		cout << " Writing grid output No = "<< myframes << " ...." << endl;
    	}
    }
    nframes++;
    return 0;
}
/*
(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
     Here, you can do whatever analysis your program requires for a frame.
	GridEps & MyEps=*MyEps0;
    Voronoi & vor=*(static_cast<Voronoi *>(data));
    MyDipole & Dip=*Dips;

	if(!symmetry){
		symmetry=new CellSymmetry(top,fr,pbc);
	}
	double ITime=0.0,t0=0.0,t1=0.0;
	if(Parameters::Input::bTiming) {
		ITime=Parallel::comm->Time();
		cout << "\n" << Parallel::comm->Get_Rank() << " MyTime is " << ITime-Mytime << endl;
		t0=Parallel::comm->Time();
	}
	symmetry->Center(nindex,cindex, natoms);
    Metric Met;
    Met(fr->box);
    Grid<1>::set(fr->box);
	Gridn<DIM>::set(fr->box);

    Dip.setCoord(Met, fr->x);
    P_t->Density(Dip);
	P_t->Filter();

	Rho_t->Density(Dip);
	(*E_t)(*Rho_t);

    if(!MyEps(*P_t, *E_t,false)) exit(1);
	if(Parameters::Input::bTiming) {
		t1=Parallel::comm->Time();
		cout << "\n" << Parallel::comm->Get_Rank() << " t1-t0 is " << t1-t0 << endl;
		cout << Parallel::comm->Get_Rank() << " --- " << fr->time << " --- --- " << fr->step << "\n";
	}
    if(Parameters::Input::freq > 0){
    	if(nframes ==0){
    		dtt=static_cast<double> (fr->time)/static_cast<double> (fr->step);
    		nstart=fr->step;
    		freq=rint(Parameters::Input::freq/dtt/Parallel::comm->Get_Size());
    		Parallel::comm->Broadcast(&freq,1);
    	}

    	if(nframes > 0) {
    		int actual=fr->step-nstart;
    		if(!(actual % freq)){
    			myframes++;
    			MyEps.setMetric(fr->box);
    			MyEps.setxRef(Dvect(0.0,0.0,0.0));
    			gmx_voroelec::fout << MyEps ;
    			cout << " Writing grid output No = "<< myframes << " ...." << endl;
    		}
    	}
    }
    nframes++;
    if(Parameters::Input::bTiming)
    	Mytime=Parallel::comm->Time();
    return 0;
}
*/

#endif /* ANALYZE_B_H_ */
